Phase separation in a binary mixture confined between symmetric parallel plates: 
Capillary condensation transition near the bulk critical point 



m 

O 
(N 



O 



I 

o 
o 



> 

p 

o 

m 



X 
J3 



Syunsuke Yabunaka^, Ryuichi Okamoto^, and Akira Onuki^ 
^Department of Physics, Kyoto University, Kyoto 606-8502 
^ Fukui Institute for Fundamental Chemistry, Kyoto University, Kyoto 606-8103 

(Dated: January 4, 2013) 

We investigate phase separation of near-critical binary mixtures between parallel symmetric walls 
in the strong adsorption regime. We take into account the renormalization effect due to the critical 
fluctuations using the recent local functional theory [J. Chem. Phys. 136,114704 (2012)]. In statics, 
a van der Waals loop is obtained in the relation between the average order parameter {tj)) in the film 
and the chemical potential when the temperature T is lower than the film critical temperature T^'' 
(in the case of an upper critical solution temperature). In dynamics, we lower T* below the capillary 
condensation line from above T^^. We calculate the subsequent time-development assuming no mass 
exchange between the film and the reservoir. In the early stage, the order parameter changes only 
in the direction perpendicular to the walls. For sufficiently deep quenching, such one-dimensional 
profiles become unstable with respect to the fiuctuations varying in the lateral directions. The 
late-stage coarsening is then accelerated by the hydrodynamic interaction. A pancake domain of 
the phase disfavored by the walls finally appears in the middle of the film. 

PACS numbers: 64.75. St,64.70.qj,68.03.Fg 



I. INTRODUCTION 

The phase behavior of fluids confined in narrow regions 
has been studied extensively [ll-[l|- It strongly depends 
on the geometry of the walls and on the molecular inter- 
actions between the fluid and the walls. Its understand- 
ing is crucial in the physics of fluids in porous media. It 
is also needed to study the dynamics of confined fluids. 

In particular, the liquid phase is usually favored by 
the walls in one-component fluids, while one component 
is pre ferentially attracted to the walls in binary mixtures 
[4|-[l0l|. In such situations, narrow regions may be filled 
with the phase favored by the walls or may hold some 
fraction of the disfavored phase. Hence, in the film ge- 
ometry, there appears a first-order phase transition be- 
tween these states, which forms a line (CCL) ending at 
a film critical point outside the bulk coexistence curve 
in the T-^oo plane fll-ls . [Tl| - [l^ . where ^oo is the reser- 
voir chemical potential I5|. We call it the capillary con- 
densation transition even for binary mixtures, though 
this name has been used for the gas-liquid phase tran- 
sition in porous media 0. Around CCL, the reservoir 
is rich in the component disfavored by the walls for bi- 
nary mixtures. With increasing the wall separation D, 
the film critical point approaches the bulk critical point. 
Crossover then occurs between two-dimensional (2d) and 
three-dimensional (3d) phase transition behaviors. 

For Ising films near the bulk criticality. Fisher and 
Nakanishi |ll| presented the scaling theory of CCL in 
the T-h plane, where h represents applied magnetic field. 
They also calculated CCL in the mean-field theory. 
Evans et al. used the density functional theory to cal- 
culate the inhomogeneous structures in pores [12l|. For 
a Lennard- Jones fiuid in cylindrical pores, Peterson et 
g^.p^ obtained steady gas-liquid two-phase patterns. 
For a lattice gas model. Binder and Landau |14| studied 



the capillary condensation transition using a Monte Carlo 
method. For a microscopic model of 2d Ising stripes, Ma- 
ciolek et al 16] found a (pseudo) CCL using a density- 
matrix renormalization-group method. For square well 
fluids in slit pores, Singh et a/.^17.] numerically examined 
the crossover from 3d to 2d. 

Recently, two of the present authors [31 calculated 
CCL near the bulk critical point using the local func- 
tional theory [l^ , which accounts for the renormal- 
ization effect due to the critical fluctuations. The low- 
ering of the film critical temperature T^^ from the bulk 
critical temperature Tc was shown to be proportional to 
(where v = 0.63) in accord with the scaling the- 
ory [ll|. Alon g C CL, our calculations jl8] and those by 
Maciolek et al [l^ showed strong enhancement of the so- 
called Casimir amplitudes [llj. Similar first-order tran- 
sitions were found between plates and colloids [11] 
in binary mixtures containing salt. 

The aim of this paper is to investigate the phase sep- 
aration in near-critical binary mixtures between paral- 
lel plates using model H and model B [l^, HI]. Here, 
phase separation takes place around CCL and the hy- 
drodynamic interaction is crucial in the late-stage phase 
separation. It is worth noting that near-critical fluids in 
porous media exhibit history-dependent frozen domains 
and activated dynamics with non-exponential relaxations 
[i^,!!^]. To gain insight into such complicated effects, we 
may start with near-critical fluids in the fllni geometry. 
Treating near-critical fluids, we may construct a univer- 
sal theory with a few materials-independent parameters, 
where D much exceeds microscopic spatial scales. 

In the literature, much attention has been paid to the 
interplay of wetting and phase separation j28l - [3l| , which 
is referred to as surface-directed phase separation. How- 
ever, simulations includin g th e hydrodynamic interaction 
have not been abundant [30, l32434| . We mention that 



Tanaka and Araki [s^ integrated the model H equations 
in the semi-infinite situation and Jaiswal et al. |34| per- 
formed molecular dynamics simulation to investigate the 
hydrodynamic flow effect between parallel plates. In our 
simulation, the order parameter ijj changes in the direc- 
tion perpendicular to the walls in the strong adsorption 
regime. Then, the dynamics is one-dimensional in an 
early stage but the fluid flow in the lateral directions ac- 
celerates the late-stage coarsenin g ev en under the no-slip 
boundary condition on the walls [2 5l [30l [33 - l33 |. 

On the other hand, Porcheron and Monson [sl] nu- 
merically studied the dynamics of extrusion and intru- 
sion of liquid mercury between a cylindrical pore and a 
reservoir. Such a process is crucial in experiments of ad- 
sorption and desorption between a porous material and 
a surrounding fluid 0. In our simulation we assume no 
mass exchange imposing the periodic boundary condition 
in the lateral directions, as in the previous simulations of 
surface-directed phase separation. 

The organization of this paper is as follows. In Sec. II, 
we will summarize the results of the local functional the- 
ory of near-critical binary mixtures in the film geome- 
try. We will newly present some results on the phase 
behavior, which will facilitate understanding the phase 
separation near CCL. In Sec. Ill, we will present our sim- 
ulation results of the phase separation with the velocity 
field (model H) and without it (model B). 



II. THEORETICAL BACKGROUND 

This section provides the theoretical background of our 
simulation on the basis of our previous paper [Tsj. The 
Boltzmann constant ks will be set equal to unity. 



with the ratio = ^o/^o being a universal number. We 
write ijj in the coexisting two phases as ±i/)cx with 



V'cx = ^cxI'tI'', 



(2.3) 



where &cx is a constant. 

We assume that the bulk free energy F including the 
gradient part is of the local functional form [T8l - l20| . 



F = j dr[/+ireC|V^p]. 



(2.4) 



In the following, we give a simple form for the free energy 
density f = fi'ip,'''). In our theory, the critical fluctua- 
tions with sizes smaller than the correlation length ^ have 
already been coarse-grained at the starting point. 



B. Coexistence-curve exterior 



Outside CX, / is of the Ginzburg-Landau form. 



(2.5) 



Here, we have omitted the free energy contribution for 
ip — 0, whose singular part is proportional to |Tp~" 
yielding the specific heat singularity. The coefhcients r 
and u in / and C in are renormalized ones in three 
dimensions. As in the linear parametric model [36(], we 
use a nonnegative parameter w representing the distance 
from the critical point in the r-ip plane to obtain 



r/r 
u/u* 
C 



(2.6) 
(2.7) 
(2.8) 



A. Ginzburg-Landau free energy 

We suppose near-critical binary mixtures with an up- 
per critical solution temperature (UCST) Tc at a con- 
stant pressure. The order parameter is proportional to 
c — Cc, where c is the composition and Cc is its critical 
value. The reduced temperature is written as 



T = (T - r,)/T„ 



(2.1) 



In our numerical analysis, the usual critical exponents 
take the following values |25|: 

a = 0.110, 13 = 0.325, 7 = 1.240, 

J/ = 0.630, 7^ = 0.0317, (5 = 4.815. (2.2) 



where Ci and u* are constants. We may set Ci = 1 

1 /2 

by rescaling ip — >■ ip without loss of generality. 
In the present case, (Ci^o)^^^V' is dimensionless. The 
constant u* is a universal number and we set u* = 
27r2/9. The fractional powers of w in Eqs.(2.6)-(2.8) arise 
from the renormalization of the critical fluctuations with 
wavenumbers larger than the inverse correlation length 
We determine w as a function of r and ip by 



w 



(2.9) 



which is equivalent to w'^ = {r + iuip"^)^^/ Ci. Thus, 
w = T for "ip = Q and r > 0, while \ip\ oc for r = 0. 

The derivative n = df /dip at fixed r denotes the chem- 
ical potential difference between the two components [lOl 
[l^,|T^, but it will be simply called the chemical potential. 
In terms of the ratio S = t/w, it reads 



At the critical composition with t > 0, the correlation 
length is written as ^ = S,ot~'^ , where is a microscopic 
length. The coexistence curve in the region r < is 
denoted by CX. The correlation length on CX is written 
as ^ = Cokl"'': where is another microscopic length 



_ 2-a + 4{l-a)S + 5aS^ 
n ~ 6[2/3 + (1 - 2/3)5]C§ 



Ciw^'ip. 



(2.10) 



On CX, we require fi = 0, which yields the equation 
2 - a -I- 4(1 - a)S' -I- 5aS'^ = for S. On CX, this gives 



S = —1/(7 or u; = —err with a = 1.714. Together with 
Eqs.(2.3) and (2.9), we obtain 

bl, = {l + <7)a^P-^/3u*C,^o. (2.11) 

W introduce the susceptibility x = x(''"iV') defined by 

TJx = = d^f/dij^. (2.12) 

For ijj ~ and t > 0, we simply obtain x(r, 0) = 
Ci^£,qT~'^. On CX, we write Xcx — %(''', V'cx)- In terms of 
the critical amplitude ratio = x(I'''|j 0)/Xcx for r < 0, 
the susceptibility on CX reads 

Xc^R^'C^'eolrr. (2.13) 

Some calculations give — 8.82 [37]. In terms of x, 
the correlation length is expressed as ^ = [CxY^^j which 
yields the critical amplitude ratio R^ — ^o/Co — 2.99 [STj . 
For r = 0, we have ^ oc |?/'|~'^^^. 

C. Coexistence-curve interior 

The interior of CX is given by \tp\ < V'cx and r < 0, 
where we need to define the free energy density / to ex- 
amine two-phase coexistence. We assume a ■i/^^-theory 
with coefficients depending only on r, where df /dip = fi 
and d^f/dip'^ — Tc/x are continuous across the coexis- 
tence curve. We then obtain 

(/ - /cx)/Tc = (V'L/8Xcx)(V'VV'cx - 1)', (2.14) 

where /cx is the free energy density on CX and Xcx is 
defined by Eq.(2.13). We also set 

C = Ccx = Cl|ar|-''^ (2.15) 

which is the value of C on CX. The renormalization effect 
inside CX is assumed to be unchanged from that on CX 
with the same r. The /i inside CX then reads 

^i/T, = (V'VV'cx - l)^/2Xcx (2.16) 

The surface tension a between coexisting bulk two 
phases is given by the standard expression, 

a = 2r,(Ccx/Xcx)'/'^L/3 

= AsTje, (2.17) 

where ^ — ^okl " is the correlation length on CX. The 
universal number As is estimated to be 0.075 in our 
model, while its reliable value is about 0.09 

D. Near-critical fluids between parallel plates 

We suppose a near-critical fluid between parallel sym- 
metric walls in the region < z < D, where D is much 
longer than any microscopic lengths. To avoid the dis- 
cussion of the edge effect, the lateral plate dimension L 
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FIG. 1: (Color online) Normalized Id profiles 'tp{z)/^pD 
(left) and ijJioc{z)D'yTc (right) vs z/D for [r /td, fioo / ^id) = 
(-2,-14.1) (top), (-20,-30.8) (middle), and (-20,-17.8) 
(bottom). Top: Adsorption-dominated profile A with ip{z) > 
0. Middle: Two profiles B and C on the capillary condensa- 
tion line with the same grand potential Q. In (b'), the two 
curves of uiodz) enclose two regions with the same area close 
to the surface tension a. Bottom: Three profiles D,E, and F 
with the same r and /loo (see Fig.3). In (c'), il increases in 
the order of D, E, and F. 



is supposed to much exceed D. The fluid is close to the 
bulk criticality and above the prewetting transition line 
[H 0, H, 0] . We use our local functional theory, neglect- 
ing the two-dimensional thermal fluctuations with sizes 
exceeding D in the xy plane. 

We scale r and i/j in units of td cx D^^^'^ and ipD oc 
jj-P/t^ ^ respectively, deflned by 

TD = {^^/Df/\ (2.18) 
= (24^/73^.*CiCo)'/V^. (2.19) 

In equilibrium theory, it is convenient to assume that 
the fluid between the walls is in contact with a large 
reservoir containing the same binary mixture, where the 
order parameter is V'oo and the chemical potential is 

Moo = m(V^oo,t). (2.20) 
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FIG. 2; (Color online) Phase diagram of a near-critical fluid 
in a film for large adsorption in the fioo/l-iD- t/td plane. 
The bulk coexistence line is given by r < and = 0. 
On its left, there appears a first-order capillary condensa- 
tion line (red bold line) ending at a film critical point, which 
is calculated from Id profiles. Displayed also are values of 
/ioo/^'D in steady two-phase coexistence in our simulation of 
a 2D X 2D x D system (x). Starting point of our simulation 
{t < 0) is also shown (o). 
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FIG. 3: (Color online) Isothermal curves in the Hoo/ij-d- 
{^)/^pD plane, which are calculated for Id profiles at t/td = 
—2, -10, and -20. For r less than its film critical value 
(= — 3.14rD), a van der Waals loop appears. Dotted parts 
of the curves for t/td = — 10 and -20 are not stable in con- 
tact with a reservoir. Points A,B,C,D,E, and F correspond to 
the curves in Fig.l. 



Here, ^oo corresponds to magnetic field h for films of 
Ising spin systems. We are interested in the case ^oo < 
(or ipoo < 0) and Vo > 0, where ipo is the value of 
ip at the walls. If equilibrium is attained in the total 
system including the reservoir, we should minimize the 
film grand potential ft. Including the surface free energy, 
we assume the form, 



n — dr wioc — Tc dS hiijj, 



(2.21) 



where the space integral J dr is within the film, the sur- 
face integral J dS is on the walls at z = and D, and hi 
is a surface field symmetrically given on the two walls. In 
addition, we neglect the surface free energy of the form 
/ dSX'^ijP' assumed in the literature [3-0] (or we consider 
the limit \ ^ oo). 

In Eq.(2.21) cjioc is the local grand potential density 
including the gradient part. 



\TcC\Vij\\ 



wioc = + T^'^cC'IV^/^r, (2.22) 
where ujs is the excess grand potential density written as 

= /(V') - /(V'oo) - Moo(V' - V^oo)- (2.23) 

Now minimization of 17 yields the bulk equation. 



l)F T 



where C — dC /dip. 
and D are given by 

V''(x,y,0) = 



(2.24) 

The boundary conditions at z = 



(2.25) 



where ^' — d'tp/dz. 

The role of hi in this paper is simply to assure 
the strong adsorption regime Vo/!'''!^ ^ (Cifo)~^^^ 
[lol [Tsl . [Toj . where ■00 is the boundary value of -0. This 
regime is eventually realized on approaching the critical- 
ity (however small hi is). In our simulation, the profile 
of ^p in the region < z < ^ is nearly one-dimensional 
depending only on z even in two phase states (see Figs. 5 
and 6). It decays slowly as (io + z)~^/'^ for < z < ^ 
[isj . where is a short microscopic length intro- 
duced by Rudnick and Jasnow Q . With the gradient free 
energy in the form of Eq.(2.22), is expressed as 

^0 ^ UCi^^y'^^i^r'^ - D{^dI^^)-'''^ (2.26) 
in terms of 'i/'o- The excess surface adsorption of in 
the region < z < is of order -0o^o ~ V'o"''^^ and is 
negligible for large 00 from P/v ^ 2, while that in the 
re gion Iq < z < ^ in oi order (^-P/'' for ^(r, 0„i) < D/2 
[sl-llOj. In the strong- adsorption regime we calculate the 
average of "0 along the z axis, 

(0) = / dz^j/D. (2.27) 



In two-phase states, depends on {x,y). 

From Eq.(2.25) it follows the estimation 
C(0'o)0'o/4- As hi/lrf^-" oo, we find 



hi ^4-''^f'iBi+B2T4'o'"' + ■■■), 

where Bi and B2 are positive constants. This is the ex- 
pression for D — > 00. In this regime, the surface free en- 
ergy in Eq.(2.21) is given by -2Bi^pI''^^A, where A = 
is the surface area. In our previous paper [18|, we ex- 
amined the film phase behavior at fixed large 0O) treat- 
ing the surface free energy as a constant. In our sim- 
ulation, we assume the boundary condition (2.25) with 
hi/C = IQllipo/D to obtain Tpo/ijjD = 14.8. 



hi 



(2.28) 



E. Capillary condensation transition 

We consider the capillary condensation transition on 
the basis of one-dimensional (Id) profiles ip = i^iz). From 
Eq.(2.25), we have the symmetry = — z). In 
the region < z < D/2, Eq.(2.24) is integrated to give 



C(V')/2 



11/2 



a;,(?A) +n 



(2.29) 



Here, IT = —A~^dVl/ dD is the osmotic pressure. It is 
the force density per unit area exerted by the fluid to the 
plates. In our case, 11 < 0, indicating attractive inter- 
wall interaction. In the Id case, it is also written as 



n = /(V^oc) - /(V^m) - Aioo(V'o 



tpm) 



(2.30) 



At the midpoint z = D/2, we set tpm = ^{.D/2). The 
fluid at the midpoint can be in the phase favored by the 
walls with V'cx due to the strong adsorption on 

the walls or in the disfavored phase with tpm = "fpoo < 0. 
Equation (2.30) indicates 11 = 2/iooV'cx in the former case 
and n = ~Tc{tpm — V'oo)^/2xcx — in the latter case, so 
n can be very different in these two cases (3^ . 

Figure 1 displays typical Id profiles of tp{z) from 
Eq.(2.29) and iviodz) in Eq.(2.22) in the range < 
z < D/2, which will be needed to explain our simu- 
lation results. Here we set (t/t/j, t/'oo/V'D; Moo /md) = 
(-2,-1.14,-14.1) (top), (-20,-1.85,-30.8) (middle), 
and (—20, —1. 81, —17.8) (bottom), where iiaa is measured 
in units of 



^iD - TJD^D oc D^/P- 



(2.31) 



Salient features in Fig. I are as follows, (i) In Fig. 1(a), tp 
is positive in the whole region with < tp > /tp£, ~ 1.20. 
(ii) In Figs. 1(b) and (b'), the fluid is on the capil- 
lary condensation line, where we give two equilibrium 
profiles B and C with the same VI. Here, B repre- 
sents an adsorption-dominated state with tpm ~ V'cx and 

< tp > /tpo = 1.99, while for C the film center is oc- 
cupied by the disfavored phase with tpm ^ —ipcx and 

< tp > /tpD = 0.187. In the right panel (b'), the integral 
of wioc(z) in the region 0<2:<D/2is the same for B 
and C. The enclosed two regions have the same area 24.7 
in units of Tc/D^ , which is close to the surface tension 

a = 29.2Tc/D'^ at this r. In addition, H a/D for B 

and H = for C [38^. (iii) In Figs.l(c) and (c'), the pa- 
rameters are those slightly below the capillary condensa- 
tion line (in Fig. 2 below). Here, there are three solutions 
with the common r and /ioo, but < tp > /tpu is 2.03 for 
(D), 0.028 for (E), and 1.20 for (F) (see Fig.3 below). If 
we perform simulation in contact with a reservoir with 
these T and /Koo, the profile D is realized at long times. 
In the very early stage of our simulation, the dynamics 
is one-dimensional and the profile F is approached after 
quenching from A in Fig. 1(a) (see Fig. 7 below). 




FIG. 4: (Color online) Phase diagrams in the r/Tn-tpm/tpD 
plane in (a) and in the t /TD-{tp) /tpn plane in (b), where 
tpm is the midpoint value. In (a) and (b), the capillary con- 
densation curves (bold lines) are calculated from Id profiles, 
where points for five t/td are those along the z axis with 
[x, y) = {D, D) and {D, 0) in the final two-phase states in our 
simulation (see Figs. 4 and 5). Bulk coexistence curve is in 
broken line. In (b) our phase separation process is illustrated 
by arrows, where the total order parameter is conserved. 



In Fig. 2, we show the capillary condensation line 
(CCL) from Id profiles located on the left of the bulk 
coexistence line in the T-fx^o plane. In our previous pa- 
per [l^ , the corresponding phase diagram was displayed 
in the r-tp^o plane. The discontinuities of the physical 
quantities across CCL increase with increasing |t| van- 
ishing at a film critical point. At this film criticality, r, 
tpao, and /ioo are calculated as 



V'oo Moo 
Td' tpo' fJ-D 



(-3.14,-1.27,-16.3), (2.32) 



where we also have {tp)/tpo = 0.989 and tp„i/ipD ~ 
—0.173. Hereafter, the chemical potential ^oo on this 
CCL will be written as ij.^^{t). Our numerically calcu- 
lated CCL is well fitted to the linear form. 



McxW/md + 16.3 = OMir/TD + 3.14). 



(2.33) 



In Fig.l(b'), the two areas enclosed by the two curves 
of u}iocD^/Tc are the same (~ aD'^/Tc). Thus, for 



TABLE I: Values of ii.oojP'D from our simulation in a finite 
2D X 2D X D system for five t /td- The corresponding values 
of on the CCL from Id profiles are also shown. 



t/td 


-8.5 


-10 


-13.3 


-16.6 


-20 


A*oo/md (finite) 


-25.8 


-27.8 


-32.2 


-36.6 


-41.4 


mS/md (Id) 


-20.6 


-21.8 


-24.9 


-27.7 


-30.8 



|t|/to 3> 1, the surface tension a and the free energy dif- 
ference per unit area are of the same order. 
See the sentences below Eq.(2.30) and the explanation of 
Fig.l(b'). For |t| ^ td, it follows the relation, 



H'cx 



-a/ijjcKD - \t/td\ 



21^-/3 



(2.34) 



Since — /3 = 0.94, the theoretical formula (2.34) is 
consistent with the numerical formula (2.33). Note that 
Eq.(2.34) is equivalent to the Kelvin equation known for 
the gas- liquid transition in pores (T], 0] . 

We have already presented a special case of three Id 
profiles in Fig.l(c) for (t/td, fi^/ f^o) = (-20,-17.8). 
In Fig. 3, we show isothermal curves in the /ioo-(V') plane, 
which are calculated from Id profiles with t/td — —2, 
— 10, and —20. The relation between /ioo and (■(/') is 
monotonic for t/td > —3.14 (above the film critical 
temperature), while it exhibits a van der Waals loop for 
t/to < —3.14 with three Id states in a window range 
A^ooi < fJ'oo < /ioo2 [Hi- Here, /iooi and /ioo2 coincide at 
the film criticality. The isothems consist of stable and 
unstable parts characterized by the sign of the film sus- 
ceptibility defined by 



Xfilm = Tc(9(?/')/5^oo)r- 



(2.35) 



In Fig, 3, points A,B,C,D,E, and F correspond to the 
curves in Fig.l. Dotted parts of the two curves of 
t/td = —10 and —20 are not stable in the presence of a 
mass current from a reservoir with common /ioo- 

Previously, some authors [H, [l3| calculated the stable 
parts of isotherms of the average density in the film versus 
the chemical potential. In our local functional theory, the 
three Id profiles can be calculated since a unique profile 
follows for any given set of r and (ip). In equilibrium 
fluctuation theory of films [ij, Xfiim is proportional to 
the variance of the order parameter fluctuations, so its 
negativity indicates thermodynamic instability. 

Furthermore, Fig. 4 gives the phase diagrams in the r- 
Tpm and T-(tp) planes. Bold lines represent the capillary 
condensation curve from Id profiles as in Fig. 2. In steady 
two-phase states in our simulation, ■(/;,„ and (ip) depend 
on X and y, so points for five r represent ipm{x,y) — 
tp{x,y,D/2) and (■(/') (x, y) with {x,y) — {D,D) and 
{D, 0) (see Figs. 5 and 6). The former line passes through 
a domain of the phase disfavored by the walls and the 
latter through the favored phase only. Phase diagrams 
similar to Fig. 4(b) have been obtained in experiments of 
the capillary condensation in porous media [2j. 
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FIG. 5: (Color online) Equilibrium two-phase state at r/ro = 
— 10 in a 2D x 2D x D system: ip/ipo (top) and ujiocD^ /Tc 
(middle) in the xy plane at z — D/2 (left) and in the yz plane 
at X = D (right). Bottom: ij/i^D (left) and wiocD^/Tc (right) 
along the z axis for {x, y) = {D, D) (blue bold line) and [D, 0) 
(red bold line), while dotted lines represent Id profiles from 
Eq.(2.29). 



III. PHASE SEPARATION DYNAMICS 

We performed simulation of phase separation in a 
Lx Lx D cell with L — 2D imposing the periodic bound- 
ary condition along the x and y axes. In this section, 
we describe phase separation realized for deep quench- 
ing. However, it was not realized for shallow quenching 
(|r|/rD < 7), for which Xfiim in Eq.(2.35) is positive. 



A. Dynamic equations and simulation method 

Supposing an incompressible fiuid binary mixture with 
a hor nog eneous temperature, we use the model H equa- 
tions [ij, [2^ . The order parameter -ip is a conserved 
variable governed by 



dtp 
'dt 



-V • {tpv) + XV 



6ip' 



(3.1) 
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FIG. 6: (Color online) Equilibrium two-phase state at r/rn = 
-20 in a 2D X 2D X D system: tp/ipD (top) and loiocD^ /Tc 
(middle) in the xy plane at z — D/2 (left) and in the yz 
plane at x — D (right). Bottom: ip/ipn (left) and ljioc /Tc 
(right) along the z axis for {x, y) — {D, D) (blue bold line) 
and (D, 0) (red bold line), which are closer to the Id profiles 
(dotted hues) from Eq.(2.29) than in Fig. 5. 



where A is the kinetic coefficient and the functional 
derivative SF/6ip may be calculated from Eq.(2.4) with 
the aid of Eqs.(2.10), (2.15), and (2.16) outside and in- 
side ex. We neglect the random source term originally 
present in critical dynamics [l^, UHl, because we treat 
the deviations much larger than the thermal fluctuations. 
The velocity field v satisfies V ■ v — and vanishes at 
z = and D. In the Stokes approximation [ioj, v is de- 
termined by 



Vpo + tpViSF/S^j) 



(3.2) 



where fj is the shear viscosity and the role of a pressure po 
is to ensure V = 0. See Appendix for the expression of 
the stress tensor in near-critical fluids and the derivation 
of Eq.(3.2). 

The kinetic coefficients A and fj should be treated as 
renormalized ones [l^, SO] (see the last sentence of 
Subsec.IIA). In the vicinity of the bulk coexistence curve, 
A may be approximated by 



where Xcx is the susceptibility on CX in Eq.(2.13) and 
is the mutual diffusion constant of the Stokes form. 



(3.4) 



with ^ = ^okl"" being the correlation length on CX. In 
our simulation, is of order tpcx at z = D/2 (see Figs. 5 
and 6 below), which supports Eqs.(3.3) and (3.4). The 
viscosity fj exhibits a very weak critical singularity and 
may be treated as a constant independent of r. 

In this paper, we also performed simulation for model 
B without the hydrodynamic interaction where ip 
obeys the diffusive equation. 



dt 5^' 



(3.5) 



The kinetic coefficient A is assumed to be given by 
Eqs.(3.3) and (3.4) as in the model H case. Then, com- 
paring the results from the two models, we can examine 
the role of the hydrodynamic interaction in phase sep- 
aration. Model B has been used to investigate surface- 



directed phase separation in binary alloys |31 1 . 

In integrating Eqs.(3.1) and (3.5), the mesh length was 
Aa; = £'/32 and the time interval width was At = 2 x 
lO^^to- The initial state was the Id profile A in Fig. 1(a) 
at t/td = —2 with small random numbers (^ 10^^) 
superimposed at the mesh points. At t = 0, we decreased 
T to a final reduced temperature. For i > 0, there was no 
mass exchange between the film and the reservoir so that 
the total order parameter J was fixed at 1.2Qtpj:,DL^ . 
We will measure time after quenching in units of 



(3.6) 



which is the mutual diffusion time in the film assumed 
to be much longer than the thermal diffusion time. We 
note that the natural time unit in bulk phase separation 
has been the order parameter relaxation time = / 
miiaiili. Here, to/h = (-D/0' = Rl\rlTD?'' - 100 
for \t/td\^ 10. 



B. Steady two-phase states 

For sufficiently deep quenching, we realized phase sep- 
aration to find a steady two-phase state at long times 
both for model H and model B. We could also calcu- 
late this final state more accurately from the following 
relaxation-type equation. 



SF /SF\' 



(3.7) 



A 



(3.3) 



where Lq is a constant and {5F /5il))t is the space average 
of SFSip at time t. Because of its simplicity, we integrated 
Eq.(3.7) with a fine mesh length of Ax — Z5/64. The data 
points in Figs. 2 and 4 and the snapshots in Figs. 5 and 6 
are those from the steady states of Eq.(3.7). 




FIG. 7: (Color online) Early-stage time-evolution of ili/tpn for 
model H (left) and model B (right) after quenching from the 
profile A in Fig.l(a) to t/td ~ —20. Shown are the profiles 
along the z axis for {x,y) = {D,D) at t/to = 0, 0.1, 0.2, 
and 0.3 (top) and along the x axis for {y,z) = (1.5D, 0.5-D) 
(bottom) at later times. In the top panels, the profile F in 
Fig. 1(c) is approached without noticeable differences between 
the two models. In the bottom panels, fiuctuations in the xy 
plane appear and coarsening is much quickened for model H 
than for model B, where t/to = 0.5, 0.6, and 0.7 for model H 
and 1.0, 1.1, and 1.2 for model B. 




FIG. 8: (Color online) Time-evolution of the normalized free 
energy decrease AF{t)/Tc for model H and model B (left) 
and v{t)to/D for model H (right) after quenching from the 
profile A in Fig. 1(a). At points (a), (b), (c), and (d) (right) 
snapshots of tp and v will be given in Fig. 9. 



In Fig. 4, the deviations of tpm and (ip) from those on 
CCL are surprisingly small, though the lateral dimension 
L is only 2D. However, in Table 1, the final tvifo-phase 
values of fioo considerably deviate from those on CCL. 
This may be ascribed to the relatively small size of the 
susceptibility Xcx for these cases. That is, if we set Xcx = 
A-^TciPd/iJ'D on CX, we obtain = 120 and 283 for 



the deviation of ^m/^D or (?/')/'(/' _d by of A^, we obtain 
that of f^oo/^J'D■ 

In Figs. 5 and 6, we display the final profiles of tp and 
wioc in the xy plane at z = D/2 (left) and in the zy 
plane at x — D (right) for t/td = —10 and —20. In 
these cases, ^ on CX is 0.078Z? and 0.051-D, respectively, 
which is of the order of the interface thickness. Displayed 
in the bottom panels are Id profiles of tp and wioc along 
the z axis for the two lateral points {x,y) = {D,D) and 
{D,0). These profiles are rather close to the Id profiles 
from Eq.(2.29) in accord with Fig. 4. 



C. Time evolution 

Both for model H and model B, early-stage time- 
evolution proceeds as follows. Just after quenching, tp 
changes only along the z axis to approach the Id profile 
at the final r with fixed {tp) (see Fig. 3). If this Id profile 
satisfies the instability condition Xfiim < 0, it follows 3d 
spinodal decomposition. On the other hand, if it satisfies 
the stability condition Xfiim > 0, it remains stationary in 
simulation without thermal noise. 

Figure 7 displays tp after quenching to t/td — —20. 
In the top panels, it is plotted along the z axis with 
(x, y) = {D, D) at t/to = 0, 0.1, 0.2, and 0.3. The velocity 
field nearly vanishes for model H, so there is almost no 
difference between the results of these models. However, 
in the bottom panels, the Id profile becomes unstable 
with respect to the fiuctuations varying in the xy plane 
for </<o ^ 0.5. The velocity field grows gradually for 
model H. In this second time range, the domain formation 
is much quicker for model H than for model B. 

In the left panel of Fig. 8, we show the free energy de- 
crease AF = Fit) — F{+Q) at t/td = —20 as a function 
oft for model H and model B. Here, Fit) is the total bulk 
free energy in Eq.(2.4) with F(+0) being its value just 
after quenching. Its decrease is accelerated with develop- 
ment of the fluctuations in the xy plane. The coarsening 
is slower for model B than for model H by about 5 times. 
In the right panel of Fig. 8, we show the characteristic 
velocity amplitude v{t) for model H, which is defined by 



v{tY = / dr\v\'/DL^ 



(3.8) 



-10 and —20, respectively. Here, if we multiply 



For t/td = -20, v{t) is equal to 0.021, 0.106, and 0.485 
at </<o ~ 0.5, 0.6, and 0.7, respectively, increasing up to 
0.8, in units of D/to. 

In Fig. 9, we show late stage snapshots of tp and v 
in the xz plane (top) and in the xy plane (middle) for 
t/td = —20 for model H, where t/to is equal to (a) 0.7, 
(b) 0.9, (c) 1.1, and (d) 1.3 after quenching. In (a) we 
can see a network-like domain of the disfavored phase. 
In (b) three domains can be seen, where the middle one 
is being absorbed into the bottom one, soon resulting in 
two domains at t/to ~ 1-0. This process gives rise to a 
dip in vit) in Fig. 8(b) since these two domains are consid- 
erably apart. In (c) and (d), furthermore, coalescence of 
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FIG. 9: (Color online) Cross-sectional velocity field v (arrows) and order parameter ^ (in gradation according the color bar) 
in phase separation for t/td ~ —20 in model H, where t/to is equal to (a) 0.7, (b) 0.9, (c) 1.1, and (d) 1.3 after quenching. 
Displayed are {vx, Vz) and t/j in the xz plane at y — D (top), those at y = 1.5D (middle), and {vx,Vy) and i/) in the xy plane at 
z = D/2 (bottom). Arrows below panels indicate typical magnitudes of the velocities, where a is the surface tension and fj is 
the shear viscosity. Final state is given in Fig. 6. 



these two domains is taking place. The arrows below the 
panels indicate the typical velocity, 0.15a/fj or 0.3a/fj, 
where a is the surface tension and fj is the shear viscos- 
ity. Note that the typical velocity in the late-stage bulk 
spinodal decomposition is given by Vc — O.la/fj psl l4l|. 
which follows from the stress balance a/R ^ QnfjVc/R 
with R ^ Vet being the typical domain length. In our 
case, from Eq.(2.17) these velocities are related as 

D/to - TjQnfiiD = {i/QTTA„D){a/fj). (3.9) 

Thus, D/to = 0.036ct/^ for t/td = -20 in Fig.9. 

IV. SUMMARY AND REMARKS 

In summary, we have examined the phase separation in 
a near-critical binary mixture between symmetric parallel 



plates in the strong adsorption regime around the capil- 
lary condensation line (CCL). Using model H and model 
B, simulation has been performed in a 2D x 2D x D cell. 
We summarize our main results. 

(i) In Sec. II, we have presented the singular free energy 
with the gradient part outside and inside the bulk coex- 
istence curve. Applying it to near-critical fluids between 
parallel plates, typical Id profiles have been given in 
Fig.l. CCL has been plotted in the t-^^o plane in Fig. 2. 
The points for the steady two-phase states from our sim- 
ulation are located in the left side of CCL. In Fig. 3, we 
have also found the van der Waals loop of isothermal 
curves in the ('0)-Moo plane, where r is smaller than the 
film critical value. The phase diagrams have been plot- 
ted in the T-^Jm and t-{i]j) planes in 3d (bulk) and 2d 
(film) in Fig. 4. The Kelvin relation (2.34) has also been 
obtained, since the osmotic pressure 11 is of order —a/D 
right below CCL [H. 



(ii) In Sec. Ill, we have first displayed tlie cross-sectional 
profiles of and wioc in steady two-phase states in Figs. 5 
and 6. The profiles along z axis for (x, y) = {D, D) 
and {D, 0) closely resemble the corresponding Id pro- 
files. For quenching to t/tjj = —20, we have examined 
time-evolution of ^. It occurs only along the z axis in the 
very early stage in the top panels of Fig. 7, where there is 
no difference between the results of model H and model 
B. Subsequently, inhomogeneities appear in the xy plane. 
The free energy decrease AF{t) = F{t) — F{+0) and the 
typical velocity amplitude v{t) defined in Eq.(3.8) have 
been plotted in Fig. 8. The velocity field considerably 
quicken the interface formation and the coarsening for 
model H than for model B. Profiles of -0 and v in the late 
stage coarsening due to the flow have been presented in 
Fig. 9, where the domain coalescence can be seen and the 
maximum velocity is of order O.lcr/fj as in bulk spinodal 
decomposition |4l|. 

We make some remarks. 

1) In the static part of our theory, we neglect the 
thermal fluctuations varying in the lateral directions 
with wavelengths longer than D. Thus this 2d transition 
exhibits mean-field behavior. In fact, the curves of V'm 
vs r and {ip) vs r are parabolic near the film criticality 
in Fig. 4. 

2) In our simulation, we soon have only two or three 
domains in the cell as in the bottom panels of Fig. 9. 
The lateral dimension L = 2D in this paper is too short 
to investigate the domain growth law in the xy plane. 
Simulation with larger L/D should be performed in 
future. 

3) From the van der Waals loop of the isothermal curves 
in the {ip)-^oo plane in Fig. 3, we may predict how phase 
separation proceeds after quenching. We have examined 
phase separation via spinodal decomposition. However, 
in real experiments, phase separation may occur via nu- 
cleation for metastable Id profiles. Note that hysteretic 
behavior has been observed in phase-separating fiuids in 
pores and has not been well explained [l|, H, [la H^l • 

4) In a number of experiments and simulations of 
surface-directed phase separation 33, 34], compo- 
sition waves along the z axis have been observed near 
the wall in the early stage. In these cases, the degree of 
adsorption has changed appreciably upon quenching. In 
the strong adsorption regime in this paper. Id dynamics 
occurs in the initial stage, but there are no composition 
waves as in the top panels of Fig. 7. 



5) The static part of this work is applicable to any 
Ising-like near-critical systems and can readily be 
generalized to n-component spin systems. In the 
dynamics, we have used model H with a homogeneous 
temperature and incompressible flows. On the other 
hand, in one-component near-critical fluids, the latent 
heat released or absorbed at the interfaces gives rise to 
signiflcant hydrodynamic flow because of the enhanced 
isobaric thermal expansion [2^ |4^. Also promising in 
future should be extension of this work to near-critical 
fiuids in porous media. 
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Appendix: Stress tensor in near-critical fluids 

In near-critical fluids, we treat slow flows with low 
Reynolds numbers. The total (reversible) stress tensor 
is given by po^ij +11^^ , where po is nearly homogeneous 
throughout the film and the reservoir. The II, ,, is the 
stress tensor due to the composition deviation j25|], 

n^ij = p^5ij + c(ViV)(Vj?/'), (Ai) 

where = d /dxi. The diagonal part p.0 is written as 

= V('5W)-/-C^|VVlV2, (A2) 

where C" = dC /dcf}. The second part in Eq.(Al) contains 
off-diagonal components relevant for curved interfaces. 
In deriving Eq.(3.2), we use the relation, 

Vj-n^^j = ipy.,{6F/5il;). (A3) 
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